For #141
I'm not quite done yet, but I wanted to ask for feedback to check that the direction is alright.

For example, the change of variables z = exp(x) looks currently like this:

eqs = [D(x) ~ α*x]
u0 = [x => 1.0]
sys = ODESystem(eqs; defaults=u0)

@variables z(t)
forward_subs  = [exp(x) => z]
backward_subs = [x => log(z)]
new_sys = changeofvariables(sys, forward_subs, backward_subs)
@test equations(new_sys)[1] == (D(z) ~ α*z*log(z))

Since there is not symbolic way (AFAIK) to get the inverse of general nonlinear functions, one has to
provide both: the forward and backward substitution (forward_subs and backward_subs).

In the test file, I also added a Riccati ODE and a 3D linear transformation as examples.

Details:

• If sys.defaults exists, then those will be transformed as well.
• The substituted variables become observables.
• Transformation can depend on time. (For transformation of initial data defaults, one needs to provide t0)
• It only works for first-order ODEs

Questions:

• Is the function interface ok? I picked this form, since it looks mathematically correct. I use expressions instead of Functions because everything else is also symbolic.

Things I will add for final PR/known bugs:

• support cases where length(forward_subs) < length(states(sys)), i.e. keep old state variables.
• add proper doc string once the interface is fixed.
• support for SDEs

Optional ideas:

• compute forward/backward substitution automatically if possible (but I don't know where to get inveres from). One could use inverse function theorem, but then the transformed ODEs are ugly/too complex.
• allow to pass known invariants, to cover cases like Riccati equations better.